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Abstract 

We investigate thermal and isothermal symmetric liquid-vapor separations via an FFT-thermal 
lattice Boltzmann (FFT-TLB) model. Structure factor, domain size, and Minkowski functionals 
are employed to characterize the density and velocity fields, as well as to understand the config- 
urations and the kinetic processes. Compared with the isothermal phase separation, the freedom 
in temperature prolongs the spinodal decomposition (SD) stage and induces different rheological 
and morphological behaviors in the thermal system. After the transient procedure, both the ther- 
mal and isothermal separations show power-law scalings in domain growth; while the exponent for 
thermal system is lower than that for isothermal system. With respect to the density field, the 
isothermal system presents more likely bicontinuous configurations with narrower interfaces, while 
the thermal system presents more likely configurations with scattered bubbles. Heat creation, con- 
duction, and lower inter facial stresses are the main reasons for the differences in thermal system. 
Different from the isothermal case, the release of latent heat causes the changing of local temper- 
ature which results in new local mechanical balance. When the Prandtl number becomes smaller, 
the system approaches thermodynamical equilibrium much more quickly. The increasing of mean 
temperature makes the interfacial stress lower in the following way: a = ao[{Tc — T)/(Tc — Tq)]^/^, 
where Tc is the critical temperature and (Tq is the interfacial stress at a reference temperature Tq, 
which is the main reason for the prolonged SD stage and the lower growth exponent in thermal case. 
Besides thermodynamics, we probe how the local viscosities influence the morphology of the phase 
separating system. We find that, for both the isothermal and thermal cases, the growth exponents 
and local flow velocities are inversely proportional to the corresponding viscosities. Compared with 
isothermal case, the local flow velocity depends not only on viscosity but also on temperature. 

PACS numbers: 47.11. Ij, 47.20.Hw, 05.70.Ln 
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I. INTRODUCTION 



Multiphase flows and heat transfers are ubiquitous in natural, industrial processes, as 
well as daily life, e.g., oil- water systems, bubble flows, petroleum processing, paper-pulping, 
and power plants, etc Therefore, establishing accurate prediction models to investigate 
the underlying physical essence of these phenomena, is of great academic significance and 
industrial practical value. However, due to the complex nature and inherent nonlinearities 
of multiphase flows, theoretical solutions are usually limited to a small class of problems 
in one-dimension and with numerous simplifying assumptions and generalizations On 
the other hand, experimental approaches for multiphase flows are generally expensive and 
some problems are still being unsolved in accurate measurement technology (e.g., interfacial 
area measurement) for this process [3']. Consequently, it is reasonable to consider numerical 
simulation, to some extent, as a primarily useful tool in studying the underlying physics of 
multiphase flows and providing some insights into understanding the kinetic process, that 
are difficult to obtain from theoretical analysis or experiments. 

Molecular dynamics (MD) is a nice microscopic approach, but it is too computationally 
expensive to access dynamic behaviors with spatiotemporal scales comparable with exper- 
iments {4]. Moreover, many macroscopic behaviors are, in fact, not sensitive to degrees of 
freedom at the molecular level. Traditional fluid dynamics does not work well for systems 
where non-equilibrium effects are pronounced, for example, multiphase system. In addi- 
tion, from the computational expenses point of view, the direct simulation of fluid behaviors 
in =uch a .y^ter. . ato a challenging work, .nee . not easy to t.aek the deto^mable 
macroscopic interfaces and the incorporate complex microscopic interactions [5,]. 

Between these two approaches, as a mesoscopic approach, the lattice Boltzmann (LB) 
method, has enjoyed substantial development and has become a very promising and versatile 
tool for simulating complex phenomena in various fields during the past two decades jo]. 



ranging from magnetohydrodynamics 



7, 



to compressible flows 



qHuJ, wave propagations 



12| . hydrodynamic instabilities 13|, llJ], etc. Apart from fields listed above, the versatile 



method is particularly promising in the area of multiphase systems [15l433j . This is mainly 
owing to its intrinsic kinetic nature, which makes the inter-particle interactions (IPI) be 
incorporated easily and flexibly and, in fact, the IPI is the underlying microscopic physical 



reason for phase separation and interfacial tension in multiphase systems. So far, many LB 
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models for multiphase flows have been proposed, among which the three well-known ones 
are the Chromodynamic model by Gunstensen et al. 15|, the pseudo-potential model by 



Shan and Chen 



16 



29( 1 , and the free energy model by Swift et al. 



■ Il7|. 



The aforementioned models have been successfully applied to study a wide variety of 



multiphase f 
wetting 20 



ow problems in science and engineering, such as contact line motion 



18 



2l|, drop breakup [22, l23|. drop collision 



phase separation and phase ordering 



15-|l7 



m |27|, 



24j |. chemically reactive fluid 



19 



32| . etc. Despite this, to date, most 



studies focus on the isothermal systems, because, in these models, only mass and momentum 
conservations are kept, hydrodynamic behaviors due to temperature fleld are not taken into 
account. However, thermal effects are signiflcant, even dominant, in many cases. Examples 
are referred to phase separations in boiling process, distillation and condensation process, and 
thermal nuclear reactor, etc. In these systems, the evolutions of the temperature fleld and 
flow fleld is spontaneously coupled with each other 
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35| . Therefore, it is a fundamental 



and essential work to develop thermal LB (TLB) models for multiphase system. But due to 
the complexity of this problem, the progress has been rather slow. 

The most obvious obstacle lies in the fact that, when the interparticle forces are incor- 
porated, how to ensure the total energy conservation becomes challenging in the discrete 
model. To overcome this difficulty, extensive efforts have been made over the past years. 
But until very recently, only a few TLB models for multiphase flows have been proposed, 
and can be roughly divided into two approaches. The flrst is the passive scalar approach 
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371]. In this approach, evolutions of the density fleld and the momentum fleld are solved 



by an isothermal LB model, while the temperature evolution is determined by an additional 
passive-scalar equation. The coupling of these two parts is through a suitably deflned body 
force in the isotherm LB equation. This approach is conceptually rather simple and as stable 
as the isothermal LB models, because the energy conservation is not explicitly implemented. 
Meanwhile, it can produce a non-ideal gas equation of state (EOS) and capture the tem- 
perature fleld. However, it should be pointed out that, in the passive scalar approach, the 
viscous dissipation and compression work done by the pressure are neglected [37 1. 

The second is the multispeed approach, which implements energy conservation by using 
larger and more isotropic sets of velocities and by including higher order velocity terms in 
;he equilibrium distribution. Examples for ideal gas include the works of Alexander et al. 



38|, Watari et al. js], Xu et al. [lo|, Q, and so on. However, applications of this approach 
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for thermal flows with high Mach number or flows with high Knudsen number still have 
some challenges. The challenges arise from the insufficient truncation in the equilibrium dis- 
tribution function and the insufficient isotropy in the discrete- velo city- mo del (DVM). In an 
alternative way, using the Hermite expansion approach, Shan et al. 39| presented a system- 
atic theoretical framework for constructing TLB models that approximate the continuum 
Boltzmann equation with higher accuracy. With the Hermite expansion approach, hydro- 
dynamic moments at various levels can be determined in a straightforward way at a given 
order of truncations of the Hermite polynomials. Almost simultaneously, similar results were 

n 

obtained by Philippi et al. [40] using a different procedure. Although the above-mentioned 
TLB models only work for ideal gas systems, they can be extended to multiphase flows by 
the extra force method. The one developed by Gonnella, Lamura, and Sofonea (GLS) 41 1 
is typical. In this model, an extra term Jfcj, accounting for inter-particle forces, is added 
into the LB equation to describe the van der Waals (VDW) fluids. From the point view of 
the IPI, it can be considered as a bottom-up approach, similar to the Shan-Chen model. To 
describe system with interfaces, gradient contributions to free energy due to the inhomo- 
geneity of fluid density are also included. Compared with the passive scalar approach, all 
observable flelds, e.g., density, velocity, temperature, and pressure are directly derived from 
the same distribution function, as in the standard kinetic theory. 

In a recent work [42], we further develop GLS model so that the total energy conservation 
can be better held and the spurious velocities can be damped to negligible scale in the 
numerical simulations. In the improved model, spatial derivatives in the convection term 
and the force term are calculated via the fast Fourier transformation (FFT) and its inverse 
(IFFT). For convenience of description, we refer to this model as FFT-TLB model. Via the 
FFT-TLB model, we study the effects of temperature and viscosity on liquid-vapor phase 
separation in two-dimensional case. It is known that, spatial domains of homogeneous phases 
evolving during spinodal decomposition (SD) show a large variety of complex spatial patterns 
and the system is globally in a nonequilibrium state. How to effectively describe and pick 
up information from such a complex system is still an open probleni. In the present work, 
besides the rheological behavior, we use the Minkowski functionals [43] to characterize the 
isothermal and thermal phase separations and conduct a comparison study on the similarities 
and differences between these two cases. 

The following part of the paper is planned as follows. The Minkowski functionals and the 
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FFT-TLB model are briefly reviewed in Sees. II and III, respectively. Simulation results 
and corresponding physical interpretations are given in Sec. IV. Sec. V presents conclusions 
and discussions. 



II. MORPHOLOGICAL CHARACTERIZATION 



In this section, we briefly review the set of statistics known as Minkowski functionals 



43|, which will be used to characterize the physical flelds in Sec. IV. Such a description 



has been well known in digital picture analysis 



the reaction-diffusion systems 



separation of complex fluids 



44| and successfully adapted to characterize 



47,^ 



shocked porous materials 



46| . and patterns in phase 



etc. 

According to a general theorem of integral geometry, all properties of a (i-dimensional 
convex set, which satisfy motion invariance and additivity, are contained in d-l- 1 numerical 
valuta. Fo. a pi.,.ed map V'(x), we consider the excursion sets of the map, deflned 
as the set of all map pixels with value of ip greater than some threshold ipth, where x is 
the position, ip can be a state variable hke density p, temperature T, or pressure P; if) 
can also be the velocity u or its components, or some speciflc stress, etc. Then the d + 1 
functionals of these excursion sets completely describe the morphological properties of the 
underlying map 'ipi^- the case of two- or three-dimensions, the Minkowski functionals 
have intuitive geometric interpretations. For a two-dimensional density map p(x), the three 
Minkowski functionals correspond geometrically to the fractional area A of the high density 
domains, the boundary length L between the the high and low density domains, and the 
Euler characteristic x- 

In this work, we probe the effects of temperature and velocity on phase separation by 
checking the density map p(x, t) and velocity map u(x, t), where time t is explicitly denoted. 
When the density p{'x,t) is beyond the threshold value pth, the grid node at position x is 
regarded as a white vertex, otherwise it is regarded as a black one. For the square lattice, 
a pixel possesses four vertices. A region with connected white (black) pixels is deflned as 
a white (black) domain. Two neighboring white and black domains present an interface or 
boundary. When the threshold contour level pth increase from the lowest density pjam to 
the highest one Pmax? the white area fraction A = NJ/N will decrease from 1 to 0, and 
the qualitative features of the patterns will vary drastically, where A^^ is the number of 
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pixels with a density larger than pth, N = x Ny is the total number of pixels, and 
A^j, are the lattice numbers along the x and y directions; the boundary length L = Nl/N 
is defined as the ratio between the pixels separating the black and white domains, and the 
total number of pixels. With the increasing of pth, boundary length L first increases from 
at Pth = Pmin, then arrives at a maximum value Lmax and, finally decreases to again at 
Pth = Pmax; the third morphological quantity is the Euler characteristic X; defined as the 
difference of the number of connected white domains and black domains N'l normalized 
hj N, X = ~ ^x^l^- -'■^ contrast to the white area A and boundary length L, the 
Euler characteristic x describes the connectivity of the domains in a purely topological way. 
It is negative (positive) if many disconnected black (white) regions dominate the image. A 
vanishing Euler characteristic indicates a highly connected structure with equal numbers of 
black and white domains. Despite having global meaning, the Euler characteristic x can 
be calculated in a local way using the additivity relation 
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481. Since the measures are 



normalized by A^, they can be used to compare systems with different sizes. 



III. FFT-TLB MULTIPHASE MODEL 



In this section, we present the FFT-TLB model for simulating thermal liquid-vapor sys 



tem. The model is a further development of the one proposed by GLS [411] • GLS introduced 
an appropriate inter-particle force term to describe the VDW fluids. Our contribution is 
to propose an appropriate FFT scheme, which is used to calculate the convection term and 
the force term. With this new model, the non-conservation problem of total energy due to 
spatiotemporal discretizations is much better solved and spurious currents in equilibrium 
interfaces are signiflcantly reduced in the numerical simulations. 

A. TLB multiphase model by GLS 

The GLS model includes the following two parts: (i) TLB model by Watari-Tsutahara 
(WT) {q]; (ii) an appropriate inter-particle force, Iki- The original WT model works only 
for ideal gas. It uses the following DVM: 

i — 1 i — 1 

vo = 0, Vfei = ffc[cos(— — vr), sin(— — tt)], = 1, 2, 3, 4; i = 1, 2. ..8, (1) 
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where subscript k indicates the k-th group of particle velocities whose speed is Vk and i 
indicates the direction of particle's speed. Different from the standard LB model, WT 
model uses a second upwind finite-difference (FD) scheme to calculate the convection term 
in the LB equation. The FD LB model breaks the combination of discretizations of space 
and time, which makes the particle speeds more flexible. The values of the speeds v^. 
may be determined in such a way that the temperature gets a large interval around the 
critical temperature Tc, under which the simulation is stable. This is of great importance for 
phase separation studies where long lasting simulations are needed to determine the growth 



behavior 



26|. 



Compared to WT model, the main contribution of GLS model is the introduction of the 
extra term I^i, which accounts for inter-particle forces 

dfki dfki 1 r „ „eq-. j /_v 

~df ' ~dv' ~ ~T " ^^'^ ' ^ ' 

where /^^^ is the local equilibrium distribution function; r is the spatial coordinate; r is the 
relaxation time related to the kinematic viscosity. The distribution function is related 
to the local density p, fluid velocity u, and temperature T through the following moments: 

pu = 5^v,,/,7, (4) 
PT = J]^(v,,-u)V.7- (5) 

ki 

Iki in Eq. ([2]) takes the following form: 

hi = -[A + BaiVkia - Ua) + (C + Cg){Vkia " Ua)'^]fll, (6) 



with 



A = -2{C + C,)T, (7) 
= ^[do,{P^ - pT) + dpKp - d^{Cd,u,)l (8) 

9 1 

+ -p^d^u^ + K[--{d^p){d^p){daUo) 
o Z 

- p{d^p){d^daUa) - {d^p){d^Uc,){dc,p)]}, (9) 
8 



C, = ^,d^[2qpT{d^T)]. (10) 

= 3pT/(3 — p) — 9p^/8 is the VDW EOS. Since the pressure is not monotonic 
in density, thermodynamic phase transition may occur in such a system. By setting 
dP'^/dp = 0, d'^P'"/dp^ = 0, we obtain the critical point = = 1. = 
Mdapdfsp - [pTd^pd^{M/T)]6ai3 - M{pV^p + |Vp|V2)<5a/3 is the contribution of density 
gradient to pressure tensor and M = K + HT allows a dependence of the surface tension 
on temperature, where K is the surface tension coefficient and if is a constant. It is worth 
pointing out that, in this model, the Prandtl number Pr = t]/k,t = ''"/2(r — g) can be changed 
by adjusting the parameter q in the term Cq. 



It has been shown that 



41| , under the Chapman- Enskog expansion, the above LB model 



recovers the following equations for VDW fluids: 

dtp + daipUa) = 0, (11) 

dt{pUa) + dp{pUaUfs + - a^p) = 0, (12) 

dtCT + dalerUa + (n^,/? - a„/3)u/3 - urdaT] = 0, (13) 

where IIq,^ = P^Sai3 + Aai3 is the non-viscous stress, and cTq/j = rj^daUji + djiUa — d^u^Sajs) + 
Cd^u^dafi is the dissipative tensor with the shear and bulk viscosities rj and C- ^toi = pT — 
9pV8 + ^ |Vpr /2 + pu^/2 is the total energy density. It should be mentioned that the force 
term also accounts for the potential energy — 9/8p^ and interfacial energy K\Vpr /2, which 
are sources of the kinetic energy. 



B. Our contribution: spatial discretization with FFT 

In this subsection, we will review our improvements to the TLB multiphase model: spatial 
derivatives in the convection term Vki'dfki/ dr and in the external force term I^i are calculated 
via the FFT scheme and its inverse. 

To illustrate the necessity, we present simulation results for a thermal phase separation 
process by various numerical schemes. Here the time derivative is calculated using the 
first-order forward Euler FD scheme. Spatial derivatives in I^i are calculated using the 
second-order central difference (2nd-CD) scheme. Spatial derivatives in convection term 
^ki'dfki/dr are calculated using the the 2nd-CD scheme, the Lax-Wendroff (LW) scheme, the 
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FIG. 1: (Color online) Variations of total energy density Aetoi{t) = etoi{t) — eto/(0) for a phase- 
separating system obtained from the GLS model with various schemes. Initial conditions are set 
as /9 = 1 + A, T = 0.85, Ux = Uy = 0, where A is a random density with an amplitude of 0.01. 
The remaining parameters are as follows: vi = 1.00, V2 = 1.90, ^3 = 2.90, = 4.30, r = 10~^, 
dx = dy = 1/256, dt = 10-^ K = 5x 10"^, H = 0, C = 0, q = -0.004. Periodical boundary 
conditions (PBC) are imposed on x and y directions. 



scheme 50|], and the fifth-order 



5l| . respectively. As a result, we 



non-oscillatory and non-free-parameter dissipation (NND _ 
weighted essentially non-oscillatory scheme (5th-WEN0) 
find that the total energy density etoi{t) is not conservative in simulations, even though it 
is in theoretical analysis (see Fig. 1). The non-conservation of energy is caused by errors of 
spatiotemporal discretizations. 

Aiming to solve the problem of energy non-conversation, we proposed a new algorithm 
based on FFT and its inverse 42|. This approach is especially powerful for periodic system 
and also provides spatial spectral information on field quantities. For completeness, let us 
start with the definition of Fourier transform of a discrete function f{xj) 



N-l 
j=0 



and its inverse 



N/2-1 



(14) 



(15) 



n=-N/2 
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In Eq. (fTSjl . k=2'Kn/L and L = NAx is the length of the system divided into equal 



segments. A general theorem of derivative based on FFT states that 52|-|54| 



f'ik) = ikx f{k), 



(16) 



where f'{k) is the Fourier transform of /'( the module of wave vector k, and i is 

an imaginary unit. The theorem provides a way to calculate the spatial derivative f'{xj), 
composed of the following steps: (i) transform f{xj) in real space into f{k) in reciprocal 
space; (ii) multiply f{k) with ik; (iii) take the inverse Fourier transform (IFT) of f'{k), 
then the spatial derivative f'{xj) can be obtained. Higher-order derivative, such as the nth 



derivative 



f{k) with (iA;)'^, 



{n > 2), can be obtained from a similar procedure only if we multiply 



(17) 



High order derivatives can be calculated from this convenient way is a main merit of FFT 
over FD schemes, otherwise, we should choose more stencils (more points) to approximate 
high order derivatives. 

The FFT approach has excellent accuracy properties, typically well beyond that of stan- 
dard discretization schemes. In principle, it gives the exact derivative with infinite order 
accuracy if the function is infinitely differentiable 53|-|56|. In our manuscript, using this 
virtue, the FFT scheme is designed to approximate the true spatial derivatives, as a result, 
to eliminate spurious velocities near the interface region and to guarantee energy conserva- 
tion. However, the trouble in proceeding in this manner is that, in many cases, it is difficult 
to ensure that the infinite differentiability condition is satisfied. For example, the function 
f'{xj) may have a discontinuity of the same character as the square wave. Then the discon- 
tinuity will induce oscillations, known as the Gibbs phenomenon. The Gibbs phenomenon 
influences the accuracy of the FFT not only in the neighborhood of the point of singular- 
ity, but also over the entire computational domain. Since the Gibbs phenomenon is related 
to the slow decay of the Fourier coefficients of the discontinuous function, it is nature to 
use smoothing pro cedures, which attenuate higher order Fourier coefficients to damp the 



oscillations 



53 



55 



57 



58|. A straightforward way is to multiply each Fourier coefficients 



by a smoothing factor 0"^, for instance, the Lanczons smoothing factor, the raised cosine 
smoothing factor, or the Fejer smoothing factor, etc 



53 



57 



58|. 
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In the recent work [42|, we presented a way to construct smoothing factors. Firstly, we 
expand k in Taylor series 

arcsin [sm{kAx / 2)] 



k 



Ax/2 

13 5 
-[sin(A;Ax/2) + - sin^(A;Ax/2) + — sin^(A;Ax/2) H sin'^(A;Ax/2) 



Aa;/2' ' ' ' 6 ' ' ' 40 ' ' ' 112 

where T{n) = t^~^e~^dt is the Gamma function, B(n) = Mod[—l + n, 2] is the Mod 
function and e{—l + n) is the unit step function. Next, in order to refrain the Gibbs 
oscillation, we should filter out more high frequency waves, or at least, damp the strengths 
of high frequency waves. Therefore, k may take the form of an appropriately truncated 
Taylor series expansion of sin(A;Ax/2). For example, k can take the following forms: 

(«) 

^^^^,^^ sin3(.A./2)/6 ^ (20) 
_ 3 3m'(tAx/2)/40 

+ aJ72 ■ 

and 

_ 5sin^(A:Ax/2)/112 
fc4-A;3 + , (22) 

and then the calculated spatial derivative is second-order, fourth-order, sixth-order, and 
eighth-order in precision, respectively. It is found that ki is consistent with the one used in 
Ref. ^|. Finally, smoothing factor for ki can be expressed as 

and the ones for A;2, ^3, and k^^ can be formulated in a similar way. 

As reported in our recent work j42], the lower-order smoothing factors, such as Oi and o"2; 
are much more effective to damp the strengthens of high frequency waves and may result 
in excessively smeared approximations, which are unfaithful representations of the truth 
physics. On the other hand, the higher-order smoothing factors, such as as and 174, can 
reserve more higher frequency waves but may not damp the Gibbs phenomenon when the 
discontinuities are strong enough, then cause numerical instability. This is especially true for 
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FIG. 2: (Color online) Variations of density Ap{t), momentum A(pu)(t), and total energy Aetoiii) 
for the phase-separating system described in Fig. 1. 

the case with shock waves and/or discontinuities. The smoothing factors should survive the 
dilemma of stability versus accuracy. In other words, they should be minimal but make the 
evolution stable. In the present study, we focus on the liquid-vapor system without shock 
waves and strong discontinuities. Therefore, the FFT scheme with higher-order smoothing 
factor, cr4, is used throughout our simulations. 

For comparisons, we verify the proposed FFT algorithm with the same problem described 
in Fig. 1 and display variations of density Ap(t), momentum A(pu)(t), and total energy 
density Aetoi{t) in Fig. 2, respectively. It is observed that, when the FFT scheme with 
(74 is adopted, variations of density and momentum nearly decrease to machine accuracy. 
For Aetoi{t), it oscillates at the beginning then goes to nearly a constant. Behaviors of 
Aetoiif) can be interpreted as follows. At the beginning of phase separation, the fluids 
spontaneously separate into small regions with higher and lower densities, and more liquid- 
vapor interfaces appear. Subsequently, spatial discretization errors induced by the interfaces 
(density gradients) arrive at their maxima, accounting for the initial oscillations. As time 
evolves further, under the action of surface tension, the total liquid-vapor interface length 
decreases owing to the coalesce of small domains, then the discretization errors, together 
with the amphtude of Aetoi{t) decrease. 
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FIG. 3: (Color online) Temperatures versus time t for the FFT-TLB scheme. 

After about 10^ time steps, the maximum derivation of etoi is only about 1.5 x 10^'', 
indicating that the FFT scheme has more advantage to guarantee energy conservation. Fur- 
thermore, we find that Aetoi{t) decreases with decreasing the initial random density A. 
When A decreases to 0.001, the maximum of Aetoi{t) will further decreases to 3 x 10^^ (not 
shown here). Numerically, this is owing to the smaller density gradients in the interface 
regions as A decreases that reduce the spatial discretization errors. Actually, A = 0.001 
(0.1% of the initial density) is enough to generate phase separation and is more appropriate. 
When A is large, or the initial temperature is far below the critical one, the initial state 
of the system is very far from the equilibrium and we may encounter large values of the 
fluid velocity in the early stage of simulations. Since the initial values of the velocity is zero 
everywhere, this process is responsible for a strong decrease of the local temperature (see 



Another interesting phenomenon is illustrated in Fig. 3. The mean and maximum temper- 
atures rise sharply at the initial period of phase separation, while the minimum temperature 
decrease significantly at first and rise rapidly at later times. The difference between the 
minimum and the maximum temperatures AT = T^ax — T^i^ arrives at its maximum at 
about t = 0.25. After that, AT decreases with time, and goes to a constant value (nearly 
vanishing) when t > 8. The reasons for behaviors of temperatures are that: at the initial 
stage, the potential energy — 9/8p^, a part of the free energy, is high, so the system will 



Fig. 3). 
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relax. During phase separation, part of the potential energy transforms into the kinetic 
energy, namely latent heat is locally released and conducts to the entire region. This is the 
main reason why temperatures are rising during simulations and the main difference from 
isothermal case, where latent heat is extracted from the system by fixing the temperature 
in all lattice nodes. Besides, viscous dissipation is another mechanism of heat generation. 
More preciously, heat is dissipated locally due to the friction between fiuid flows when the 
fluid velocity is different from zero. After phase separation, interfaces forms and the fluid 
velocities go to zero everywhere. Then the kinetic energy transforms into thermal energy 
totally. 

In our recent work, the FFT-TLB multiphase model has been validated successfully by 



two sets of typical benchmarks 42|]. Simulation results demonstrated that the FFT-TLB 
model can capture both qualitatively and quantitatively the interface properties in accord 
with the VDW theory. Besides that, with the new model, spurious velocities near the liquid- 
vapor interface are signiflcantly reduced, and, as a result, phase diagrams of the liquid- vapor 
system obtained from simulations are more consistent with that from theoretical calculations. 

IV. SIMULATION RESULTS, RHEOLOGICAL AND MORPHOLOGICAL 
CHARACTERIZATIONS 

When a system is suddenly quenched into the two-phase region, the original single phase 
becomes unstable, then phase separation occurs through the formation and the subsequent 
growth of domains. Eventually, the system arrives at a new equilibrium state. In the 



past few decades, this phenomenon has been extensively studied 
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11, 



21, 



34 



35 



60 



65| . by theoretical derivations, experiments, and numerical simulations. Among others, the 
most signiflcant flnding is the domain growth law, which states that, at late times, the 
characteristic domain size R{t) grows as a power with time t, R{t) ~ t". The value of 
exponent a is believed to be universal, depending only on the growth mechanism, and has 
been well known in isothermal system, a = 1/2 and 2/3 for high and lower viscosities. 



respectively 271, l6ll, |62|, |65| . However, behaviors of phase separation with temperature fleld 
are far from clear. The aim of this section is to clarify effects of temperature dynamics on 
both the rheological and morphological behaviors of phase separation. 
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A. Patterns for isothermal and thermal cases 



Simulations for isothermal and thermal phase separations are performed on lattices with 
X Ny = 512 X 512 nodes. PBC are imposed on both directions. Here, we only consider 
symmetric mixtures, namely we set liquid:vapor mass fractions to 1:1, for which at late times 
these domains will form a bicontinuous structure with sharp interfaces 66|, |67|. Therefore, 
the initial conditions are set as follows: 

(p, u^, Uy, T) = (1.042 + A, 0.0, 0.0, 0.9), (24) 

where 1.042 is the mean density of liquid and vapor at T = 0.9, and A is a random density 
noise with an amplitude of 0.001. Parameters are set to be r = 10^^, At = 10"^, K = 
5 X 10^^, Ax = Ay = 1/256, and others are unchanged. Density distribution patterns at 
representative times t = 0.4, 1.0, 2.5, and 8.0 are shown in Fig. 4 for isothermal case (see 
Figs. 4I(a)-l(d)) and thermal case (see Figs. 4II(a)-II(d)). For the isothermal case, after 
about 25000 time steps, the fluid has begun to separate spontaneously into small regions 
with higher and lower densities. As time evolves, the small domains merge with each other 
and larger domains appear under the action of surface tension at t = 0.4. From patterns at 
t = 0.4, 1.0, and 2.5, as excepted, higher and lower densities domains evolve in an equal way, 
leading to an interwoven bicontinuous pattern. The growth of domains continues at t = 8.0, 
and, eventually, the system will reach a completely separated state for a large enough time. 

Compared with configurations in the isothermal case, several distinctive differences can 
be found in the thermal case: (i) the average size of domains in each case tends to increase in 
an effort to decrease the interfacial energy, while at the same moment, in the isothermal case, 
it is bigger than its counterpart, which demonstrates that the domains grow faster in this 
case; (ii) for the isothermal case, interfaces between vapor and liquid are much clearer, which 
shows that the interfaces in this case are much narrower; (iii) density difference between the 
maximum and minimum densities Ap = Pmax — Pmin in the isothermal case is much larger 
than the one in thermal case, indicating that phase separation in this case is deeper; (iv) 
contrary to interpenetrating bicontinuous structures formed in the isothermal case, isolated 
and nearly circle vapor droplets suspending in the liquid phase are appeared in thermal 
case. These differences are interesting and meaningful. In the following subsections, we will 
analysis these differences with the help of rheological description and Minkowski functionals. 
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FIG. 4: (Color online) Snapshots in two processes of phase separations. The temperature is fixed 
at T = 0.9 in process I. See Figs. I(a)-I(d). The initial temperature in process II is T = 0.9. See 
Figs. II(a)-II(d). The relaxation time is fixed at r = 10"'^ in the two processes. The time t = 0.4, 
1.0, 2.5, and 8.0 in (a), (b), (c), and (d), respectively. The lattice size here is 512 x 512. 
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FIG. 5: (Color online) Spherically averaged structure factor S{k,t) versus wave number k for the 
procedures shown in Fig. 4. Figure (a) is for the isothermal case and Fig. (b) is for the thermal 
case. In each figure, S{k,t) at times t = 0.25, 0.3, 0.5, 0.6, and 0.8 are shown in the inset, while 
those at times t = 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, and 8.0 are shown in the main frame. 

B. Rheological characterization 



In order to further quantify the results shown in Fig. 4, time evolution of the circularly 
averaged structure factor S{k,t) is employed, which is defined as the Fourier transform of 
the density-density correlation function. For a discrete system, it can be stated as 



S{k,t) 



J][p(x,t)-p(t)]e^'^- 



/N, 



(25) 



where k = {2tt /N){mi + nj) is the wave vector in the reciprocal space with m = 1,2, ...A^x, 
n = 1,2, ...Ny. S{\<.,t) is further smoothed by averaging over an entire shell in k space to 
obtain the circularly averaged structure factor 

S{k,t) = J2S{k,t)/J2^. (26) 

k k 

In Fig. 5, we present the time evolutions of S{k, t) for isothermal case in (a) and thermal 
case in (b), respectively. All curves in Fig. 5 can be roughly divided into two different time 
regimes: the SD stage and the domain growth (DG) stage. From Fig. 5(a), at early times, 
such as t = 0.25 and 0.3, we observe that the peak in S{k, t) increases in height without the 
position of the peak changing in time. This behavior is indicative of the initial sharpening 
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of domains, without detectable phase separation taking place. In the second stage, the peak 
of S{k,t) increases in height and shifts to smaller wave number, indicating the coarsening 
of domains. At i = 0.5, 0.6, and 0.8, we observe the appearance of a second peak in S{k,t), 
which merges with the main peak later on. This behavior manifests that there is more than 
one typical domain size at that moment. From t = 6.0 onwards, the peak seems to stop 
drifting to the left but only oscillates in amplitude, which means that the finite size effects 
are pronounced. 

Similar results can also be found in the thermal case. Nevertheless, careful comparisons 
of these two cases will show you some distinctions: (i) the first stage continues up to t = 0.8, 
which is longer than the isothermal case. The existence of temperature field significantly 
decelerate the speed of domain formation, an effect which has also been seen in Fig. 4; 
(ii) over the period from t — 4.0 to t — 6.0, the peak of S{k, t) only varies in height but 
very little in wave number. This phenomenon is usually observed at the initial stage of 
phase separation, leading us to think that the system has steered to a new SD stage before 
reaching the finial late time stage. Essentially, during this stage, the dynamics is mainly 
making the interfaces thinner while the average domain sizes barely change; (iii) at the same 
time, the peak of S{k, t) in isothermal case is much larger than the one in thermal case, but 
the corresponding wave number is much smaller, which demonstrate that both the density 
difference between the two phases and the characteristic domain size are much larger in the 
isothermal case. These results agree with the information obtained from Fig. 4. 

Next, the characteristic domain size R{t) is used to further describe the kinetic process 
quantitatively. R{t) is derived from the inverse first moment of S{k,t), 



In Fig. 6, we show the growth of R{t) versus iterations for, r = 10~^ in (a), and r = 10~^ 
in (b), in a log- log scale. In each figure, the top and bottom scatter symbols correspond to 
the simulation results for the isothermal and thermal cases, respectively. Straight lines in 
each plot are linear fits of the simulation results. Discarding both the early time transient 
regime and the very late time regime, where finite-size effects are pronounced, we find, for 
isothermal case, the behaviors of R{t) during the DG stage are R{t) ~ for r = 10~^, 
and R{t) ~ for r = 10""^. These results are in good agreement with the generally 
accepted theoretical predictions of R{t) ~ i^/^ and R{t) ~ t"^^^ at high and low viscosities 




(27) 



k 



k 
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FIG. 6: (Color online) Domain growths in isothermal and thermal systems with r = 10~^ in (a) 
and T = 10~^ in (b). The squares and diamonds are for results from FFT-TLB simulations. Lines 
are shown in each plot to guide the eyes. 



by the Allen-Cahn theory 68|] and LB models 27|, |6l|, l62|, |65|]. But for the thermal case, the 



growth exponents decrease to 0.46 for r = 10~^, and 0.58 for r = lO""^, respectively. This 
can be regarded as another proof for our conclusion, which states that domains grow faster 
at lower temperature. In addition to the above differences, another piece of information also 
deserves our attention. For the isothermal case with r = 10~^, the SD stage lasts about 
for 25000 time steps (see the upper horizontal solid line in Fig. 6(a)). Nevertheless, for 
the thermal case, it lasts for 80000 time steps (see the lower horizontal dash dot line in 
Fig. 6(b)). Similar results can also be found in the case with r = 10^^. These findings 
suggest that, compared to the isothermal case, the SD stage is significantly prolonged by the 
existence of temperature field. In the following part, the morphological functionals are used 
to verify similarities and differences between the two cases, and the corresponding physical 
interpretations are given. 



C. Morphological characterization and physical interpretations 

1. Similarities and differences 

To perform Minkowski functionals analysis for the density map, we choose a density 
threshold pth and pixelize the map into high density regions (with p > pth) and low density 
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FIG. 7: (Color online) Minkowski measures for the procedures shown in Fig. 4. The left column 
is for the isothermal case and the right column is for the thermal case. 

regions (with p < pth)- Figure 7 shows the time evolutions of Minkowski measures for the 
procedures shown in Fig. 4. From Fig. 71(a) for the isothermal case, we see that, when 
Pth = 0.40, the white area fraction A keeps nearly 1.0 during the whole procedure shown 
here, which means no local density is lower than 0.40 in the system up to t = 8.0. However, 
when the threshold increases to 1.70, the white area fraction A keeps nearly zero during 
the whole process. Thus, no local density is higher than 1.70 in the system. As a direct 
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consequence of SD, A increases with time t when pth > 1-0, while it decreases when pth < 1-0. 
Consequently, most curves (except for the uppermost curve for pth — 0.40 and the lowermost 
curve for pth = 1.70) toward the horizontal central line from about t = 0.25. Afterwards, 
at the DG stage, the curve for pth = 1-05 (mean density of the system) overlaps with the 
horizontal central line and other curves are symmetric to it. The outer two curves are for the 
cases with pth — 0.45 (the red ball) and pth = 1.60 (the dark yellow hexagon), respectively. 
We mention that these two densities are just the equilibrium densities of vapor and liquid 
at T — 0.9. Obviously, the outer two curves mark the reach of the correct equilibrium 
state. As time evolves, proportion between the two curves decreases, as a consequence of 
the coarsening of the high/low density domains. Moreover, the proportion between any two 
curves can be conveniently obtained from Fig. 71(a). 

Now we go to the second and the third Minkowski measures, the boundary length L and 
Euler characteristic x- As shown in Fig. 71(b), for each case, L increases sharply to its 
maximum at about t — 0.25, then decreases slowly. The first increase and the subsequent 
decrease in L are due to the appearance of the liquid- vapor interface during the SD stage 
and the following decrease in interfacial area during the DG stage, respectively. At the SD 
stage, when pth < 1.0, x decreases to be evidently less than zero, which indicates that the 
number of domains with p < pth increases. On the contrary, when pth > 1.0, % increases 
to be evidently larger than zero, which indicates that the number of domains with p > pth 
increases. These results show that the phase separation process is in progress. We mention 
that, at about t — 0.25, the case with ptu — 0.45 has the minimum Euler characteristic and 
the case with pth = 1.60 has the maximum one, but the two cases get the minimum boundary 
length L. These results indicate the following information: for the first case, many scattered 
black domains with pt/t < 0.45 appear in the high density background with pth > 0.45, while 
for the second case, the high density domains with p^^ > 1-60 are scattered in the low density 
background with pth < 1.60. These domains are so small that the total boundary length is 
nearly zero. From Figs. 4I(c)-I(d), we observe that the density maps show highly connected 
structures with nearly equal and very small numbers of black and white domains. Hence, 
the Euler characteristic x keeps close to zero in the DG stage (see Fig. 71(c)). 

From Figs. 711(a)-ll(c), for thermal case, one can also distinguish two different stages. At 
early times {t < 0.8), due to the growth of density fluctuations and the build up of interfaces, 
the density area fraction A belonging to the liquid phase increases, while the one belonging 
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to the vapor phase decreases. The changes also result in the increase in boundary length 
L (see Fig. 711(b)). The appearance of liquid- vapor interfaces has an additional effect. 
They separate the system with disconnected minority domains. As a result, the absolute 
value of Euler characteristic x increases in the SD stage. In contrast to the first stage, as a 
direct consequence of the coalescence of relatively small domains, the characteristic length 
scale increases but the number of domains decreases, therefore, the DG stage (t > 0.8) is 
characterized by the decrease in L and x- 

In the end of the second stage, an interesting phenomenon occurs. There are two small 
proportions for p < 0.65 and p > 1.35 during the second stage and reach their maxima at 
about t = 3.0 (see Fig. 711(a)), but are gradually diminishing afterwards. This phenomenon 
shows that a recombination process is taking place owing to the increasing temperature 
that interrupts the original process and forces the system evolves to a new equilibrium state 
decided by the variable temperature. 

Compared with figures shown in Fig. 7, main differences between these two cases are 
analysed and listed as follows: (i) for isothermal case, the domain with a density between 
[0.45, 1.60] only accounts for 20% at t = 4, and decreases further with time. But for thermal 
case, the domain with a density between [0.75, 1.25] reaches to 50% of the whole domain 
at t = 8. The difference demonstrates that the separation depth in thermal case is much 
shallower, while the interface width is much wider, which can be clearly seen in Fig. 711(b); 
(ii) the maxima of L and x can be used to mark the transition from the SD stage to DG 



stage 



481]. The transition time for isothermal case is about 0.25, but for thermal case, it 



increases to about 0.8. The result further confirms our conclusion that: phase separation 
occurs faster in the isothermal case. From another point of view, this conclusion can also be 
obtained from the A{t) curve. For most cases, after the initial quick changing period, the 
changing of A with time t shows a slowing down. The slope of the A{t) curve corresponds 
approximately to the speed of phase separation. For the same density threshold, the slopes 
of A{t) curves in the two cases are quite different. For example, when pth = 0.65, the A{t) 
curve decreases sharply in isothermal case, while for the thermal case, it decreases much 
more slowly. So we can say that, in isothermal case, the process of phase separation is much 
faster; (iii) connectivity of patterns in isothermal case is much better than that in thermal 
case. This feature can be achieved from the evolution of x- X decreases enormously and 
almost vanishes at about t = 0.5 in the isothermal case, but it is negative for the vapor 
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FIG. 8: (Color online) Portion of the density and temperature gradient distributions at i = 0.8 in 
(a) and t = 7.0 in (b) for the thermal case with r = 10^^. The lengths of the temperature gradient 
vectors are magnified by 400 times in (a) and 4000 times in (b). 



structure until about t = 6.0 in the thermal case, which is consistent with density patterns 
in Fig. 4 (see Fig. 411(d)). 



2. Physical interpretations of the prolonged SD stage and the lower growth exponent in thermal 
case: effects of temperature and viscosity 

In Sec. IV, we find, compared to isothermal case, the SD stage is significantly prolonged 
and the growth exponent is lowered in thermal case. In this subsection, effects of temperature 
and viscosity are investigated to provide proper interpretations. 

Firstly, in Fig. 8, we display the density and the corresponding temperature gradient 
distributions at two representative times for the thermal case with r = 10^'^. To illustrate the 
structure of the temperature gradient fields clearly, the lengths of the vectors are multiplied 
by 400 in (a) and 4000 in (b), respectively. As shown in Fig. 8(a), many tiny droplets and 
bubbles appear in the system, and the temperature gradient vectors are toward the droplets. 
Thus, the local temperatures within droplets are slightly higher than the mean temperature 
of the system, while the local temperatures within bubbles are slightly lower than the mean 
temperature. With the separating process, the local temperatures in the two phases deviate 
more from the mean temperature and an overshoot phenomenon is observed. This procedure 
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continues up to an extent, after which the local phases with high (low) temperatures partly 
begin to transform back from liquid (vapor) to vapor (liquid). In this way, both the local 
high temperatures and low temperatures approach the mean temperature, and the system 
approaches thermodynamical equilibrium quickly at lower Pr number (Pr = 0.1 for r = 
10^^). This process is evident by Fig. 8(b), where there is no determinate relationship 
between temperature gradients and liquid (vapor) domains. Moreover, the temperature 
difference between the highest and lowest decreases to 0.014. The system approaches to 
thermodynamical equilibrium so quickly that the temperature difference becomes so small 
during the phase separation process. Therefore, in this case, temperature can not be regarded 
as an ideal physical quantity to describe this process. 

In another way, we employ enthalpy and latent heat to describe this process, and the 
enthalpy is defined by 

h = e + P/p, (28) 

where e = pT — 9/8p^ + KlVpf /2 is the internal energy density including the gradient 
contribution. The difference of enthalpy between two states determines the latent heat 

Lh = h2- hi. (29) 

Figure 9 shows the density and enthalpy distributions at t = 6.0 for the thermal case with 
r = 10^'^. Comparison of the two figures illustrates that the enthalpy of vapor is relatively 
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FIG. 10: (Color online) Distribution of latent heat (a) = ht=s.o ~ ^t=3.o ^-nd the corresponding 
distribution of density change (b) Ap = pt=%.o — Pt=3.o for the thermal case with r = 10^'^. 



higher than that of hquid. In order to study the dynamic characteristics of the pattern, we 
show the spatial distribution of latent heat = /it=8.o ~ ^t=3.o in Fig- 10(a) and density 
change Ap = pt=s.o — Pt=3.o between these two states in Fig. 10(b). Careful observations 
between these two figures suggest that when the latent heat Lh is positive, the corresponding 
density difference Ap is negative. Droplets (Bubbles) absorb latent heat and evaporation 
occur simultaneously. Subsequently, the density decreases. A negative Lh corresponds to 
an increase of density, then the droplets (bubbles) have a coagulation trend. It should be 
noted that, owing to the transformation of potential energy into thermal energy, the total 
latent heat is released during the whole process. The released heat conducts over the entire 
region rapidly at low Pr number, and increases the mean temperature of the system (see 
Fig. 11). While in isothermal system, latent heat is extracted from this system by fixing 
the temperature in all lattice nodes. 

As well, another piece of information can be obtained from Fig. 11. The mean tempera- 
ture scarcely grows at the first stage due to no detectable phase separation is taking place 
and no remarkable latent heat is released. Subsequently, in the next stage, the temperature 
rapidly increases to 0.97 and, later, keeps almost zero growth. Afterwards, under the almost 
unchanged temperature, phase separation evolves in accord with the isothermal case. 

So far, we have not discussed in detail the surface tension. For a planar interface, it can 
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FIG. 11: (Color online) Time evolutions of the mean temperature and released latent heat for the 
procedures shown in Figs. 4 II(a)-II(d). 
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FIG. 12: (Color online) Time evolution of surface tension for the procedures shown in Fig. 4. 
Figure (a) is for the isothermal case and Fig. (b) is for the thermal case. 



be computed from the following formula 69l-l7l| : 



or from the VDW theory 72474] 



where 



$* = p*^ - p*T*[ln{l/p* - 1) + 1] - p 



*2 



(30) 



(31) 



(32) 
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e = T* ln(l/p: - 1) - p:T7(1 - pD + 2p:, s = v,l 



(33) 



p* = pb, T* = bT/a, with a = 9/8 and 6 = 1/3 in this model. It is pointed out that Eq. (I3T!) 
is especially convenient, since it can be evaluated directly without determining the density 
profile. Now, we calculate surface tension with Eq. (I30l) for both the isothermal and thermal 
cases from profiles along the x-axis and, at the same time, calculate theoretical values from 
Eq. fl3Tl) . These results are plotted in Fig. 12. It is shown that, in the isothermal case, after 
the formation of liquid-vapor interface, the surface tension keeps nearly a slightly oscillating 
constant around the exact value. While in thermal case, the surface tension is much smaller 
than the theoretical one before interfaces are well formed (before the mean temperature 
reach to 0.96). After that, it decreases obviously with the increase of temperature, and can 
be verified in the following form: 



where ao.g is the surface tension at T = 0.9. The increasing temperature lowers the density 
gradient, as well as the surface tension, which is the driving force for diffusive growth. As a 
result, domains grow more slowly than in the isothermal case. 

Essentially, during the whole process, compared to the isothermal case, two competition 
mechanisms exist. The first one is heat generation and conduction mechanism, or temper- 
ature rising mechanism. The release of latent heat results in a rising temperature and the 
rising temperature results in a new dependence of pressure-density. In other words, it leads 
to a new local mechanical balance. The second one is the hydrodynamic flow generation 
and development mechanism, or liquid-vapor equilibrium mechanism, decided by viscosity, 
diffusivity of the fluid, etc. They compete and influence with each other, deciding the finical 
morphology jointly. 

Besides the temperature effects, we now consider how the hydrodynamic flows influence 
both the morphology of the phase separating system and the growth exponent. Figure 13 
shows the time evolution of the high velocity |u| area fraction for the procedures shown in 
Fig. 4. From it, also, two stages can be found, corresponding to the nearly-zero value of 
white area fraction for all |u|^^, the rapidly increase and the subsequent slowly decrease. 
For isothermal case, |u| G [0,0.09], while for thermal case, |u| G [0,0.02]. This implies that 
velocities are not only damped by viscosity but also by the rising temperature. The maximum 
velocity in the isothermal case can reach to 0.09 or even higher. Consequently, in contrast 



a = ao.9[(T,-T)/(T,-0.9)]3/^ 



(34) 
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FIG. 13: (Color online) High flow velocity area fraction A versus time t for the procedures shown 
in Fig. 4. Figure (a) is for the isothermal process and Fig. (b) is for the thermal case. 
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FIG. 14: (Color online) High flow velocity area fraction A versus time t for two phase separation 
procedures. Figure (a) is for the isothermal process and Fig. (b) is for the thermal case. Here, the 
relaxation time is set to be r = 10~^. 

to thermal case, hydrodynamic effects can not be totally neglected in the isothermal case. 
The appearance of larger flow velocities offers more opportunities for coalescence between 
domains. Under the action of diffusion and hydrodynamic flows, a faster DG process is taking 
place, and a bigger growth exponent can be observed. Figure 14 shows the same trend of the 
hydrodynamic flows, when r decreases to 10~^. Due to the lower viscosity, velocities are more 
sufficiently developed. Figures 13 and 14 demonstrate that, in the thermal case, compared 
to the thermodynamic and diffusion mechanisms, hydrodynamic flows are less important 
than that in the isothermal case and, therefore, can not be regarded as a dominant factor 
governing the growth exponent. 
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FIG. 15: (Color online) Time evolution of each part of total energy for the procedures shown in 
Figs. 4II(a)-II(d). 

In Fig. 15, we display the time evolution of each part of total energy for the procedures 
shown in Figs. 4II(a)-II(d), which presents a clear image about energy evolution during 
phase separation. It is found that, pT and — 9p^/8 are the main parts of total energy, 
and evolve in the opposite way. Kinetic energy and surface energy are much smaller than 
the two former ones. The maximum of K\V /2 corresponds to the appearance of nuclei 
and formation of small domains. Afterwards, it decreases gradually due to the increasing 
temperature and the decreasing interfacial area. The macroscopic kinetic energy pu^ j1 is so 
small that the viscous dissipation induced by it can be neglected. Therefore, in thermal case, 
compared to latent heat, the effects of kinetic energy on temperature are less important. 



V. CONCLUSIONS AND DISCUSSIONS 



Thermal and isothermal symmetric liquid-vapor separations are simulated via the FFT- 
TLB method. Structure factor, domain size, and Minkowski functionals are used to describe 
the density and velocity fields and, at the same time, to understand the configurations and 
the kinetic processes. Simulations and physical analysis present the following scenario for 
the thermal phase separation. When the separation starts, many tiny droplets and bubbles 
appear in the system. The local temperatures within droplets are slightly higher than the 
one within bubbles. With separating, neighboring droplets (bubbles) coalesce and the mean 
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domain size increases. The local temperatures in the two phases deviate more from the 
mean temperature. This procedure continues up to a stage, after which the local phases 
with high (low) temperatures partly begin to transform back from liquid (vapor) to vapor 
(liquid). In this way, both the local high temperatures and low temperatures approach the 
mean temperature, and the system approaches thermodynamical equilibrium. 

Simulation results also indicate that the phase separation in thermal and isothermal cases 
can be generally divided into two stages: the SD stage and the DG stage. Different from 
the isothermal case, the SD stage is significantly prolonged, and different rheological and 
morphological behaviors are induced by the variable temperature field in the thermal case. 
After the transient procedure, both the thermal and isothermal separations show power-law 
scalings in the domain growth; while the exponent for thermal system is lower than that 
for isothermal system. With respect to the density field, the isothermal system presents 
more likely bicontinuous configurations with narrower interfaces, while, the thermal system 
presents more likely configurations with scattered bubbles. 

Compared with the isothermal case, heat creation, conduction, and lower interfacial 
stresses are the main reasons for the differences in thermal system. Latent heat, is released 
during the separating process, which is the main reason for the rising temperature. The 
changing of local temperature results in new local mechanical balance. When the Prandtl 
number becomes smaller, the system approaches thermodynamical equilibrium more quickly. 
The increasing local temperature has an additional effect. It makes the interfacial stress 
lower. This behavior in simulations is quantitatively verified by the theoretical formula, 
cr = (To [(Tc — T)/{Tc — To)]^/^, where Tc is the critical temperature and uo is the interfacial 
stress at a reference temperature Tq. Besides thermodynamics, we find that the local viscosi- 
ties also influence the morphology of the phase separating system. For both the isothermal 
and thermal cases, growth exponents and local flow velocities are inversely proportional to 
the corresponding viscosities. Compared with isothermal case, the local flow velocities in 
thermal case not only depend on viscosity but also temperature. In future studies, we will 
increase the depth of separation which the FFT-TLB model can undergo, and investigate 
quantitatively how the Prandtl number affects the separation procedure. 
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